#include <stdio.h>
#include <stdlib.h>
#include "mkl_lapack.h"
#include "cublas_v2.h"
#include "mkl_lapack.h"
#include "dept150_common.h"
#include "cuda.h"
#include "cuda_runtime_api.h"


#define BLOCK_SIZE 16
#define PRINT_ALL1
int LPPSADB_cu(int n, double* a, double* b, double *E, double *El, double *cond, double ea, double eb);
int LPPSADB_cu(int block_size, int n, double* a, double* b, double *E, double *El, double *cond, double ea, double eb);

int LPPSADB_cu(int n, double* a, double* b, double *E, double *El, double *cond, double ea, double eb)
{
	LPPSADB_cu(BLOCK_SIZE, n, a, b, E, El, cond, ea, eb);
}

int LPPSADB_cu(int block_size, int n, double* a, double* b, double *E, double *El, double *cond, double ea, double eb)
{
	double *dev_matrix;
	double *host_matrix_block;
	int pivot_row;
	int mkl_info;
	cublasHandle_t handle;
	cublasStatus_t stat;
	cudaStream_t stream;

	host_matrix_block = (double*)malloc(block_size*block_size*sizeof(double));

	gpuDeviceInit(0);
	//checkCudaErrors(cudaStreamCreate(&stream));
	
	stat = cublasCreate(&handle);
	if(stat != CUBLAS_STATUS_SUCCESS)
	{
		printf("CUBLAS INITIALIZATION FAILED");
		return EXIT_FAILURE;
	}

	checkCudaErrors(cudaMalloc((void **)&dev_matrix, n*n*sizeof(double)));
	
	stat = cublasSetMatrix(n, n, sizeof(double), a, n, dev_matrix, n);
	if(stat != CUBLAS_STATUS_SUCCESS)
	{
		printf("CUBLAS UPLOAD MATRIX FAILED");
		return EXIT_FAILURE;
	}
#ifdef PRINT_ALL
	printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");
	cublasGetMatrix(n, n, sizeof(double), dev_matrix, n, a, n);
	print_matrix(n, n, a, n);
#endif
	pivot_row = 0;
	mkl_info = 0;
	double alpha_1 = 1;
	double alpha_m1 = -1;
	while(pivot_row + block_size < n)
	{
		//printf("ITER %d\n", pivot_row);
		//copy pivot block to host and factorize it using mkl
		stat = cublasGetMatrix(block_size, block_size, sizeof(double), &dev_matrix[IDF2(pivot_row, pivot_row,n)], n, host_matrix_block, block_size);
		if (stat != CUBLAS_STATUS_SUCCESS)
		{
			printf("CUBLAS DOWNLOAD MATRIX FAILED");
			return EXIT_FAILURE;
		}
		
		dpotrf("L", &block_size, host_matrix_block, &block_size, &mkl_info);
		if(mkl_info)
		{
			printf("BLOCK FACTORIZATION FAILED %d", mkl_info);
			return EXIT_FAILURE;
		}
		
		cublasSetMatrix(block_size, block_size, sizeof(double), host_matrix_block, block_size, 
			&dev_matrix[IDF2(pivot_row, pivot_row,n)], n);
		
		//find A21
		stat = cublasDtrsm_v2(handle, CUBLAS_SIDE_RIGHT, CUBLAS_FILL_MODE_LOWER, CUBLAS_OP_T, CUBLAS_DIAG_NON_UNIT, n - pivot_row - block_size,
			block_size, &alpha_1, &dev_matrix[IDF2(pivot_row, pivot_row, n)], n, &dev_matrix[IDF2(pivot_row + block_size, pivot_row, n)], n);
		if(stat != CUBLAS_STATUS_SUCCESS)
		{
			printf("CUBLA DTRSM ERROR %d %d", stat, (n - pivot_row - block_size));
			return EXIT_FAILURE;
		}

		//update A22
		stat = cublasDsyrk_v2(handle, CUBLAS_FILL_MODE_LOWER, CUBLAS_OP_N, n - pivot_row - block_size, block_size, &alpha_m1, 
			 &dev_matrix[IDF2(pivot_row + block_size, pivot_row, n)], n, &alpha_1, &dev_matrix[IDF2(pivot_row + block_size, pivot_row + block_size, n)], n); 
		if (stat != CUBLAS_STATUS_SUCCESS)
		{
			printf("CUBLAS DSYRK FAILED");
			return EXIT_FAILURE;
		}

		pivot_row+= block_size;
	}
	//copy last block to host and factorize it
	stat = cublasGetMatrix(n - pivot_row, n - pivot_row, sizeof(double), &dev_matrix[pivot_row* n + pivot_row], n, host_matrix_block, block_size);
	if (stat != CUBLAS_STATUS_SUCCESS)
	{
		printf("CUBLAS DOWNLOAD MATRIX FAILED");
		return EXIT_FAILURE;
	}
	int cur_block = n - pivot_row;
	dpotrf("L", &cur_block, host_matrix_block, &block_size, &mkl_info);
	if(mkl_info)
	{
		printf("BLOCK FACTORIZATION FAILED %d", mkl_info);
		return EXIT_FAILURE;
	}
	
	//copy matrix to host
	cublasGetMatrix(n, n, sizeof(double), dev_matrix, n, a, n);
	//todo set the rest of the matrix
	for(int i =0; i < cur_block; i++)
	{
		for(int j =0; j < cur_block; j++)
		{
			a[IDF2(pivot_row + i, pivot_row + j, n)] = host_matrix_block[IDF2(i,j, block_size)];
		}
	}
}
